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We present the revised "Meudon" model of Photon Dominated Region 
(PDR code), presently available on the web under the Gnu Public Licence at: 
http: / / aristote.obspm.fr/MIS. General organisation of the code is described down 
to a level that should allow most observers to use it as an interpretation tool with 
minimal help from our part. Two grids of models, one for low excitation diffuse 
■ clouds and one for dense highly illuminated clouds, are discussed, and some new 

results on PDR modelisation highlighted. 
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Subject headings: Astrochemistry - ISM: general - ISM: molecules - Methods: 
numerical. 



1. Introduction 

The recognition of the importance of photodissociation processes and their relation to 
atomic to molecular transitions in interstellar clouds started with Bates & Spitzer (1951). 
However, the first numerical quantitative model computing the various atomic and molecular 
abundances was presented by Black & Dalgarno (1977). These authors focused their study 
on the diffuse cloud toward the bright star ( Oph described as a "standard" diffuse cloud 
for which extensive observational information was available, including the first detection of 
molecular hydrogen in various rotational levels derived from absorption transitions with the 
Copernicus satellite (Morton 1975). Further impetus to study such environments has been 
provided on the one hand by the detection of infrared emission of molecular hydrogen, ionised 
and neutral carbon, neutral oxygen and the so called aromatic bands at 3.3, 7.6 and 11 /xm. 
The availability of the facility cryogenic grating spectrometer with the 91cm telescope of the 
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Kuiper Airborne Observatory has allowed to detect for the first time, the far infrared emission 
lines of [01] (63, 146 //m), [CII] (158 /im), [Sill] (35 jum), and excited rotational lines of CO in 
luminous galaxies, planetary and reflection nebulae. These observations have emphasized the 
importance of Photon Dominated Regions and were interpreted in the context of theoretical 
models such as those of Tielens & Hollenbach (1985a,b) and Le Bourlot et al. (1993a). The 
ISO (Infrared Space Observatory) observed bright nebulae such as in Orion (Abergel et al. 
2003). 

These almost neutral astrophysical regions where ultraviolet photons penetrate and af- 
fect the physical and chemical properties are described now as Photon Dominated Regions 
(PDRs). These include hot gas close to HII regions, diffuse and translucent clouds (Galactic 
and inter-galactic), and the envelopes of dark clouds where star formation takes place. They 
have been extensively described in the review article of Hollenbach & Tielens (1999). Addi- 
tional impulse came after the launch of the FUSE (Far Ultraviolet Spectroscopic Explorer) 
satellite in 1999 which can detect H 2 and HD in absorption in a variety of Galactic and 
extragalactic environments (Rachford et al. 2002; Tumlinson et al. 2002; Bluhm et al. 2003; 
Boisse et al. 2005; Lacour et al. 2005). The perspective of Spitzer and Herschel reinforces 
the need to describe as accurately as possible the physical conditions of photon-dominated 
environments and clarify the possible diagnostics. Various groups have developed PDR codes 
at different stages of sophistication and a workshop has been held in the Lorentz Center in 
Leiden in March 2004 to make detailed comparisons between the different numerical codes. 
Roellig et al. (2005) describe the results obtained for different test cases. 

The present paper delineates the recent implementations performed in our PDR model 
(Le Bourlot et al. 1993a), labelled as the "Meudon PDR code", which is available at the 
address http://aristote.obspm.fr/MIS. The aim is to gather in a single paper the physical 
questions and the numerical answers which have been implemented so far and to provide to 
the observers a numerical tool to interpret their observations. We also want to emphasise 
remaining issues to consider, such as the inclusion of data and/or processes which are poorly 
known. The organisation of the paper is the following: Sects. 2 to 7 describe the physics 
and give technical details on the organisation of the code (2: general features, 3: grains, 
4: radiative transfer, 5: excitation, 6: thermal balance, 7: chemistry). Some representative 
examples are found in Sect. 8 and 9, and Sect. 10 is our conclusion. Various appendices 
highlight more specific points. 
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2. General features 
2.1. Overview 

The model considers a stationary plane-parallel slab of gas and dust illuminated by an 
ultraviolet radiation field (described in Appendix C) coming from one or both sides of the 
cloud (the two intensities can be different). It solves, in an iterative way, at each point 
in the cloud, the radiative transfer in the UV taking into account the absorption in the 
continuum by dust and in discrete transitions of H and H 2 . Explicit treatment is performed 
for C and S photoionisation, H 2 and HD photodissociation and CO (and its isotopomeres) 
predissociating lines. The model also computes the thermal balance taking into account 
heating processes such as the photoelectric effect on dust, chemistry, cosmic rays, etc. and 
cooling resulting from infrared and millimeter emission of the abundant ions, atoms and/or 
molecules. Chemistry is solved, and abundances of each species are computed at each point. 
The excitation state of a few important species is then computed. The program is then able 
to calculate column densities and emissivities/intensities. 

One major restriction is the strict steady-state approximation, so that model results 
cannot be compared directly to observations of rapidly evolving regions. However, the time 
scales of photoprocesses are modest at lower extinction and/or high radiation fields and 
shorter than the typical two body chemical reaction time scale in diffuse environments. The 
time scale given by the H 2 photodissociation is typically 1000/x years at the edge of a cloud 
(with x the FUV radiation strength - see Sect 4.2). The steady state approximation is then 
satisfactory. 

2.2. Variables and parameters 

As modellers, we can define the parameters which describe the system and which can 
be tuned at our will. As a consequence, some variables will adjust themselves under those 
choices and the constraints of physical laws. 

Our first hypothesis is that each cell of gas is small enough for all physical quantities 
to be constant in it, but large enough for statistical mean to be meaningful. We can thus 
speak of quantities such as "kinetic temperature" (Tk) as a function of position. Note that 
this single very general hypothesis rules out some interesting problems, such as the presence 
of shocks. 

The two most important physical quantities considered are density and temperature. 
Here density is the total number of hydrogen nuclei in all forms % = n(H) + 2 n(H 2 ) +n(H + ) 



-4- 



in cm -3 . To keep enough flexibility, they may be either parameters or variables depending on 
the user's choice. In the first case, a complete cloud structure can be specified (by analytic 
expressions, or in a "profile" input file). Temperature becomes a variable if thermal balance 
equations are solved, and density itself becomes a variable if some kind of state equation is 
used. The most usual cases are to solve for thermal balance and use either a constant density 
or a constant pressure case. Other variables have a better defined status and are described 
in Table 1. 

Model parameters are much more numerous. They fall into two main categories: 

Astrophysical parameters define the cloud characteristics and environment and are dis- 
played in Table 2. They are tunable over several orders of magnitude as will be discussed 
below. 

(Micro)physical properties, on the other hand, should have well defined values (e.g. 
chemical reaction rate coefficients) within the range of considered variables but are in fact 
also parameters and displayed in Table 3. Indeed, some properties suffer from different 
uncertainties and offer an opportunity of choice to the modelist. These uncertainties have 
to be remembered and taken into account to assess the sensitivity of results. 

2.3. Structure of the code 

Running a model consists of three distinct steps. First, choices have to be decided on 
a set of parameters appropriate to the goal of the computation. This is done by filling the 
input data files, which allows handling very different conditions. It is generally meaningless 
to try and use some published graph in the hope of interpreting some specific observational 
data. A much better way consists in running and adapting a model to its particular needs. 
Second, the physical structure of the cloud is computed from the input file parameters and 
the result is saved as a data base in a binary file. This is the computationally intensive part. 
Finally a post-processor code allows the user to dig into this binary file to extract quantities 
of interest and perform the physical and chemical interpretation. It is only during this final 
step that integration of local properties along the line of sight is achieved and "observational" 
quantities (such as, e.g., column densities, line intensities) are computed. 

During the second stage, and to maintain reasonable computation times, many aspects 
of "real" clouds have to be idealised. While making changes to the code over the years, 
some approximations have been refined and more accurate processes included, but some 
very basic assumptions can not be overcome without writing a new code from scratch. The 
most fundamental of these are that this a one-dimensional, plane parallel, steady-state code: 



- 5 - 



• All quantities, either input parameters or computed variables depend at most on one 
spatial coordinate (some are even constant but are probably varying in a true cloud, 
as e.g. grain characteristics). 

• The spatial coordinate is taken as the dust optical depth in the visible r v computed 
perpendicular to a semi-infinite plane. Code design makes it nearly impossible to turn 
it e.g. to a spherical geometry. 

• In all physical equations, terms involving are set to 0. This implies that we compute 
the state of the cloud after infinite time has passed. In practical, infinity means longer 
than any characteristic time of the problem. This holds for all physical processes (e.g. 
formation and destruction times of molecules), and for the external environment of the 
cloud. 

These conditions place limits on the kinds of objects that can be modelised. For example, 
it is not possible to compute intensities from a strictly edge on PDR since the line of sight 
in that direction is infinite. It is also not very wise to compare emissivities to those from a 
very young pre-stellar object, since time dependent effects are bound to dominate. 

The precise structure of the modelled cloud is described in Appendix A. 

Coupling of various physical processes (such as radiation transfer and chemistry, thermal 
and chemical balance) requires using various numerical methods during the computations. 
As the radiation field decreases as the optical depth ry increases, the energy density inside 
the cloud is reduced accordingly. Therefore, physical conditions at a given point depend 
more on parts of the cloud closer to the nearest edge than on more shielded parts. This 
property is used to compute physical quantities from one edge up to deeper parts. However, 
it is not possible to reach a complete self consistent solution in one step. First, if radiation 
is allowed to come from both sides, the computation eventually reaches regions where the 
influence of not yet computed parts of the cloud is far from negligible. Second, when the 
gas opacity is taken into account in the UV radiative transfer, it is mandatory to know 
in advance the level populations, which depend on the radiation field to be computed. So 
the solution is reached through an iterative process, where an approximate structure of the 
whole cloud is saved at step % and used to compute a refined one at step i + Convergence 
properties are shown in Appendix B. 

At each iteration, physical processes are treated in turn in the following order: 
1. Radiative transfer in the UV 



2. Chemistry 
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3. Level populations and thermal balance 

As radiative transfer depends on the populations of the various species and chemistry is 
sensitive to temperature, a second level of iteration is needed at each point into the cloud. 
Chemical equations themselves are solved using a Newton-Raphson scheme which is itself an 
iterative method. Thus, there is a hierarchy of embedded iterative processes, each of which 
has to proceed until a convergence criterium is reached. 

Finally, physical quantities are supposed to vary (almost) continuously with ry, the 
visual optical depth. However, these variations may be very steep (e.g. at the H/H 2 tran- 
sition), and the optical depths where gradients of the abundances are high, are not known 
in advance. Also, the intrinsically exponential character of radiation attenuation with depth 
implies that very different scales must be taken into account. We deal with those constrains 
by using an adaptative spatial grid with logarithmic variations. The first point is at r v = 
(with perfect vacuum behind in place of ionised gas), and the second at t v = 1CT 8 . r v 
grows from that value by factors of 10 until some variations of the main physical quantities 
are obtained. From this point, the spatial step is kept short enough to ensure that relative 
variations of all important quantities are kept below a predefined threshold. 

Although unphysical, discontinuities may appear at some particular value of the optical 
depth r v . This comes from the possible occurrence of multiple solutions to the steady state 
equations for some set of parameters (Le Bourlot et al. 1993b). As r v increases, physical 
conditions vary smoothly until the point where one of the solutions vanishes. If the model 
result follows this branch of solution, it will "jump" to the other branch with a discontinuity. 
In a "real" cloud, diffusion or turbulent mixing would smooth those discontinuities. This 
point is further discussed in Sect. 7. 

3. Grain properties 

Grain properties are involved in, at least, three important physical aspects : 

• They determine the cloud extinction curve, which is needed for UV radiative transfer. 

• They may catalyse some chemical reactions, and are accountable for all H 2 formation 
in standard galactic conditions. 

• They take part in thermal balance through photo-electric heating, and weak collisional 
coupling with the gas. 



-7- 



In an ideal model, these three contributions should derive from a uniquely defined distribu- 
tion. However, this is impossible to achieve within the present knowledge of grain physical 
properties and more empirical approaches are required. We have extracted three sets of 
parameters which are mainly involved in the radiative UV transfer and the thermal balance 



- the grain size distribution, which is a power law function, following the analysis of 
Mathis et al. (1977) and referred as the MRN law: 

dn(a) oc a a da 

where dn(a) the number of spherical grains per unit of volume to have a radius between 
a and a + da and a = —3.5. 

- the UV extinction curve as a function of the wavelength. A standard galactic expres- 
sion can be used as well as specific functions depending on the environment. We use the 
parametrization of Fitzpatrick & Massa (1990): 

E X -y x 2 

c 1 + c 2 x + c 3 — -o— — + c 4 F{x) 



Eb-v ° (x 2 — ul) 2 + x 2 ^ 2 

with F(x) = 0.5392 (x - 5.9) 2 + 0.05644(x - 5.9) 3 
if x > 5.9 /im" 1 and elsewhere, with x — 1/A in jimT 1 . 

The six parameters ci,C2, C3, C4, 7 and yo are taken to be mainly a function of Ry = 
A y /E(B-V). 

- grain absorption cross sections taken from Draine & Lee (1984) and Laor & Draine 
(1993). These values are determined from experimental analysis of extinction properties of 
spherical particles of graphite and silicates for a sample of radii between 1 nm and 10 /im. 
Values corresponding to radii smaller than 50 nm are considered to mimic PAH (Poly Aro- 
matic Hydrocarbons) properties. 

Currently, graphite and silicates are mixed in constant proportion within the code. 
Additional grain properties such as the dust to gas mass ratio, the volumic mass and the 
mean distance between adsorption sites on a grain etc. which affect the reactivity on dust, 
are displayed in Table 4. 

The physical effects resulting from the gas-grain interaction have been described in Le 
Bourlot et al. (1995b) and the model follows the corresponding prescription. 

The present choice of parameters differs somewhat from other conventions found in other 
PDR models. It is indeed often assumed that a single parameter, the effective UV continuum 
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absorption cross section (a = ira 2 Q cxt (l00nm) n g /n}i) at 100 nm allows to take care of the 
continuum UV flux attenuation. Sternberg & Dalgarno (1989) introduced an effective grain 
absorption cross section of 1.9 x 10~ 21 cm 2 , supposed to hold in the far-ultraviolet spectrum 
where H 2 absorption occurs. It is interesting to note that, with our choice of parameters, 
the mean geometrical cross section of grains is equal to 1.1 x 10 -21 cm 2 with the Mathis 
size distribution of grains and 5.8 x 10~ 22 cm 2 with a unique value of the radius of 0.1 /im. 

4. U.V. Radiative transfer 
4.1. Strategy 

The resolution of the complete radiative transfer equation is an exceedingly time con- 
suming computation because : 

• The whole spectral range from the far ultraviolet (FUV) to millimetre wavelength 
needs to be considered, implying a variety of processes to consider 

• The spectral resolution in discrete lines requires a high number of grid points (even 
using a variable grid size) 

• radiative processes are coupled to chemical and thermal processes in a non-linear way. 

Therefore the problem has to be split, and various approximations introduced to render it 
tractable. Our main assumptions are the following: 

• Incoming UV radiation is decoupled from outgoing IR and millimetre radiation. 

• There is no internal UV source function, so that transfer at UV wavelengths is scat- 
tering or pure absorption followed by emission. Photodissociation and photoionisation 
resulting from secondary UV photons generated by cosmic rays on the molecular gas 
are not explicitly calculated. Scaling factors to the cosmic ray ionisation are introduced 
for each species with the appropriate dissociative or ionised channel as computed by 
Gredel et al. (1989). 

• Redistribution of radiation effects are neglected except for the anisotropic UV scat- 
tering by grains. In particular, cooling lines are treated within the "on-the-spot" 
approximation: Photons are either re-absorbed where they are emitted or escape from 
the cloud. Therefore, level populations can be computed from detailed balance using 
only local quantities. 
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• During the first stage of the computation (building the data base) the source function 
in various millimeter or submillimeter cooling transitions is obtained with an escape 
probability approximation (de Jong et al. 1980) for computing the line opacity; how- 
ever, no approximation is done in the integration over optical depth at line center for 
each line. 

• During the post-processing phase, when temperature, abundances and populations are 
determined at each point inside the cloud, the full transfer equation can be solved 
along a line of sight to obtain line profiles and integrated intensities (see Fig. (1)). 

• We do not solve the thermal balance of the grains resulting from the IR continuum 
emission of the dust, nor the emission from PAH. 



4.2. UV Continuum 

The external UV radiation field is described in Appendix C. By default, the Draine 
(1978) radiation field is used and is scaled by a multiplicative factor x given in the entry 
file. Its attenuation into the cloud is computed by solving the radiative transfer equation on 
a variable grid of wavelengths using the spherical harmonics method (Flannery et al. 1980; 
Roberge 1983). 

This equation reads: 

where \x = cosd with 9 the angle between the direction perpendicular to the slab of gas 
and the direction of propagation of the beam of light (see Fig. (1)), s is the curvilinear 
abscissa in the direction perpendicular to the slab of gas. I{r,fi) is the specific intensity 
in ergs -1 cm -2 sr" 1 A -1 at the position r in the direction given by /i where k,\ is the total 
absorption coefficient (gas plus dust) in cm -1 , and <x A the dust scattering coefficient (cm -1 ). 
p(fx, jj,') is the angular redistribution function by dust of the radiation field. Gas extinction 
in selected lines of H and H 2 may be included by using as an independent variable the true 
opacity: 

dr x = -ds (k% + k,x + <7 A ) 

where is the dust extinction and is the gas absorption coefficient, all in cm -1 and at 
the wavelength A. 
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The contribution from dust is easily computed from (see also Appendix A) 



Ry E~Q_y 



and 

C 1 
ds = 1.086 — — dry 
Ry Tin 

Extinction curves ( can be selected on a per object basis from an easily expandable 
data base, using the analytical fits of Fitzpatrick & Massa (1986, 1988, 1990). Many authors 
have given the extinction curves using this same set of 6 parameters. The Ry factor can be 
determined observationally (Cardelli 1994; Patriarchi et al. 2001, 2003) and has a mean value 
close to 3.1. If unknown, in diffuse gas, the standard mean Galactic parameters reported in 
Fitzpatrick & Massa (1990) are used. 



4.3. UV discrete absorption 

Discrete absorption in the UV also occurs and is mainly due to atomic H Lyman lines 
and molecular H2 electronic transitions within the Lyman and Werner system bands. Pho- 
todissociation of H 2 takes place via these discrete transitions coming from specific rotational 
levels (J=0, 1, etc.). These lines may become very wide due to saturation effects. The set of 
transitions to introduce in the radiative transfer can be chosen from none to any at the user's 
choice and is determined from the value of the rotational quantum number of the lower level 
of the UV transition. A huge computation time can result if many lines are included. Adding 
more and more lines does not lead to a simple scaling between the number of lines and the 
computing time. CPU time is proportional to the number of points in the wavelength grid 
and line overlap tends to limit its growth as the number of lines increases. 

For the lines taken into account in the radiative transfer the method is the following. 
For a given absorption line from lower level I to upper level u, the contribution to opacity 
at wavelength A is computed from 

where Az/d = ^ \j + "^turb * s the nne Doppler width, Tk the gas kinetic temperature, 
fturb the turbulent velocity, rii the density of the lower level population, f iu the oscillator 
strength of the transition I — > u. 
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The Voigt profile is defined by 

a f + °° e- y2 dy 



*J-oo « 2 + (e-y) 2 



where a = 4?r ^ , T is the upper level inverse life time and £ = The Voigt function 

is computed using Wells (1999) 1 . Various tests have shown that it is faster to compute the 
contribution of all lines at all wavelengths than to try and devise position and wavelength 
dependent tests to include only significant contributions. 

The transfer equation then becomes: 



lM?^- = I(r,iM)-S(r,iM) (1) 

with: 

S(r, 1*) = ^ f ' pfo //)/(r, //) d» (2) 

Differences with the method of Roberge (1983) come from the dependence of the effective 
albedo uo with r. This dependence results from the contribution of the gas to the UV 
absorption. Details are given in Appendix D. 

We use a variable grid in wavelengths to maintain accuracy while limiting memory and 
CPU requirements. The grid is built during the initialisation phase in two stages: 



• first, all "required" wavelengths are listed. They include various values, either phys- 
ical (ionisation thresholds, line positions, etc.) or imposed by input files (grids from 
external files, etc.). 

• Then the grid is completed by choosing points in the profiles of each line (if any) using 
increasing steps from each centre of line. A test integral is evaluated, and the grid 
refined until adequate convergence is reached. 



Various numerical tools have been designed to manipulate the radiation field, allowing dy- 
namical memory allocation. All parameters are tunable, but not included in the user input 
file: this part is considered too technical for promotion to general "user space" . 

As an example, a portion of the radiation field inside a small clump (Ay = 0.5 mag) of 
density 100 cm~ 3 and irradiated by the standard ISRF on both faces is presented in Figure 



Source code and corrections available at http://www.atm.ox.ac.uk/user/wells/voigt.html 
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(2). Note that line transfer is a CPU time consuming process. Including no lines in the 
radiative transfer, only H 2 lines from J=0 and 1, or H 2 lines from J=0 to J=3, CPU times 
are in a ratio 1 : 10 : 50. 



4.4. Other lines 

The H 2 lines above the threshold mentioned in the previous paragraph and the CO 
predissociating lines (as well as the ones of their isotopes since D, 13 C or 18 are introduced in 
the model) are not included in the complete radiative transfer described above. We compute 
self-shielding effects by using the approximation of Federman et al. (1979) (hereafter FGK 
approximation). It gives a fair approximation of the optical depth at the centre of the 
line and of the resulting self-shielding from a knowledge of the abundance and temperature 
along the line of sight, known from a previous iteration. Note that absorption is computed 
using the local (ry dependent) radiation field, which includes absorption by H 2 strong lines. 
Protection of CO and HD by H 2 is then taken into account. 

4.5. UV escape probability 

From Flannery et al. (1980) we have an elegant mean to compute the escape probability 
of UV photons produced within the cloud, e.g. within H 2 cascades. Those photons are 
emitted isotropically so that the fraction which reaches one side of the cloud is a transfer 
problem similar to the one described in the previous section. The cloud is split into two 
parts at the position f of the emitting molecule. Boundary conditions (i.e. impinging 
specific intensities) are set to 1 at f and to at r = and at r = r max (where r max is 
the maximal optical depth of the cloud), and Eq (Dl) is solved for the two half-clouds. 
The probability for a photon emitted at f to reach (respectively r max ) is then the ratio 
J(0)/J(f) (respectively J {r max ) / J {f)) from that sub-system. Those escape probabilities are 
computed once per global iteration, thus reducing the CPU requirements. 

The main drawback of this method is the "on the spot" approximation: photons that 
do not escape the cloud are assumed to be absorbed where they were emitted. This underes- 
timates the energy loss from that point and changes slightly the level populations. However, 
lifting that approximation would require a full treatment of radiative transfer (including the 
line source function) which is out of the scope of our model. 
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5. Excitation 

Level populations are computed for a number of species (see Table 5). Detailed balance 
is solved explicitly including collisional transitions, radiative decay, interaction with the 
cosmological background, and state to state formation and destruction processes. 

We follow de Jong et al. (1980) to evaluate a resonant photon escape probability. Solving 
the non-linear coupled steady state equations provides both level populations and cooling 
rates through radiative decay. Dust extinction is neglected for infrared and millimeter lines. 

Note that for one-sided models, the escape probability towards the interior of the cloud 
is strictly zero. As we use the "on the spot" approximation, this results in a modification of 
the local source function that may modify significantly the populations. Thus, it is always 
better to compute two-sided models if an estimation of the real size is known, even if the 
radiation field is very different on both sides. See Sect. 6.2.1 for a quantitative example. 

After all iterations have converged, the whole cloud state is known, and the radiative 
transfer equation can be solved to compute a line emissivity in a given direction and the line 
profile. Fig. (3) shows the profile of C + fine structure line at 158 /xm for the same cloud as in 
Fig. (2). The profile is essentially Gaussian, dominated by turbulent velocity (here 2kms~ 1 ), 
and with very little saturation (optical depth at line centre: 0.24). The effect of radiation 
attenuation by H 2 lines is barely detectable here: the integrated line intensity decreases from 
3.8 10~ 6 to 3.3 10~ 6 ergcm~ 2 s -1 sr -1 . This results from the reduced ionisation of C inside 
the cloud due to shielding by H 2 lines (C column density varies from 2.1 10 14 to 2.9 10 14 
cm -2 ). 

6. Thermal Balance 
6.1. Heating 

Heating results from impinging UV photons and cosmic rays, but several different pro- 
cesses are involved to convert either directly or indirectly the energy input into kinetic energy. 

6.1.1. Photo-electric effect on grains 

We follow the development of Bakes & Tielens (1994), but add some significant upgrades. 
Grain absorption coefficients are taken from Laor & Draine (1993). Integration over radiation 
field is performed using the local field as computed from the radiative transfer (see Sect. 4). 
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The charge distribution is computed for each grain size. Finally, integration over size is 
performed. 

6.1.2. H 2 formation on grains 

H 2 formation releases about 4.5 eV. However that energy is unevenly distributed be- 
tween grain excitation, internal energy and kinetic energy of released molecules. The branch- 
ing ratios are not known and probably depend on conditions in the cloud and the nature of 
the grain itself. Our treatment is based on an equipartition argument, but with technicalities, 
and H 2 formation as discussed further in Le Bourlot et al. (1995c). 

• First, the internal energy distribution of newly formed H 2 on grain can be chosen 
among several options. This results in some internal energy E- mt which is computed 
before each run. Note that some hypotheses may lead to most available energy to be 
in E int , precluding any heating. 

• Secondly, the fraction of energy remaining in the grain is set to one third of H 2 forma- 
tion energy at most, provided there is enough available. 

• The remaining fraction is set to kinetic energy. 

Using standard options, this results in effectively 1.5 eV released in the gas by molecule 
formation. But the code adjusts itself automatically to other prescriptions, providing one is 
available. 

The initial population of H 2 upon formation is the same throughout the modeled cloud. 
Our standard option is to distribute an energy of Ei nt = 1.478 eV with a Boltzmann distri- 
bution. Note that this is not equivalent to using a Boltzmann factor of Tg = E [nt /k. Given 
n ma x, the highest included level of H 2 , we have: 

YZ=i En 9n exp (- J?;) 

E- m t — 7 \ 

ET = T 9n exp (- 

where E n is the energy of level n and g n its statistical weight. Including all levels of H 2 , 
we find that T B = 8734 K is required to recover E- mt /k = 17322 K (1.478 eV), while using 
T B = 17322 K leads to E int /k = 28237 K. This may induce very significant differences in 
the computed gas temperature in regions where active formation/destruction of H 2 has a 
significant contribution to the thermal balance. 
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The initial ortho to para ratio upon formation is set to 3. That value may be changed 
if needed, which has a large influence on the final ortho to para ratio in regions where 
photodestruction occurs on time scales smaller than the conversion ones. However, there is 
no firm experimental indication yet that this occurs on interstellar grains. 

6.1.3. Photodissociation and Photoionisation 

The kinetic energy of the atomic and/or molecular fragments released in photodissoci- 
ation and photoionization may also contribute to the heating of the gas. As these processes 
may also lead to internal excitation which is rarely known, we have introduced different as- 
sumptions, depending on the species involved. The kinetic energy released in the photodisso- 
ciation of H2 has been computed by Abgrall et al. (2000) for each discrete level involved. The 
corresponding mean kinetic energy is about 1 eV, in agreement with the previous estimation 
of Stephens & Dalgarno (1973). We have assumed the same energy release for HD, CO and 
its isotopes. The energy input (erg cm -3 s -1 ) is then obtained by multiplying this energy 
by the photodissociation probability (s" 1 ) and the abundances of the species involved. For 
other species, we have assumed a mean photon energy of 13 eV. In atomic photoionization 
(C, S, Fe, ..) we assume that all the available energy is transferred in kinetic form. When 
molecules are produced, we assume that one third of the available energy is in kinetic form 
and participates to the heating. Abundant species such as H 2 , C, CO give the main con- 
tribution but this term is never dominant. Similar assumtions are used for secondary UV 
photons. 

6.1.4- Cosmic rays destruction 

The heating rate due to cosmic ray ionisation is also poorly known. Black (1987) esti- 
mated that this contribution should be between 4 and 6 eV. We assume 4 eV per ionization 
by cosmic rays. 

6.1.5. Chemical equilibrium 

We calculate the energy balance of chemical reactions from the variation in formation 
enthalpies between the products and the reactants. This term contributes to the heating 
of the gas and becomes important in dark and cold regions when UV photons are almost 
completely absorbed. We note that the estimated contribution is probably an upper limit : 
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indeed, the main fraction is due to dissociative recombination reactions which involve large 
exothermicities and exothermic ion-neutral reactions. We assume that molecular species are 
produced at thermal equilibrium and do not consider branching ratios in excited states. The 
actual source of this heating is cosmic rays which produce the ions. 

6.1.6. Other heating sources 

Additional heating terms coming from macroscopic motions, such as turbulence or shock 
waves may be introduced via empirical formulae such as those derived by Black (1987). 
Although we have done various tests, these are not included in the standard version of the 
code as most of the physics is still poorly known or can not be expressed via simple analytical 
formulae. 



6.2. Cooling 

Cooling comes mainly from discrete radiative transitions in lines of various species fol- 
lowing collisional excitation. 

6.2.1. Fine- structure transitions 

These are described in Sect. 5 and Table 5. Note that cooling is underestimated for one- 
sided models since the "on the spot" approximation forbids photons from escaping towards 
the "rear" side of the cloud. 

As an example, we show in Figure 4 the temperature computed at the edge of a cloud 
of density n H = 10 3 cur 3 illuminated by a radiation field of x — 10 3 . The total depth of the 
cloud varies from Ay = 1 to Ay = 100 mag, and we use either radiation coming from one 
side only or from both sides (with x — 1 on the far side in this case). Three features are 
worth noticing: 

1. For Ay > 5 mag, j£ dge is constant, and 1 side and 2 sides models agree. The small 
difference (a few degrees) can be ascribed to numerical convergence effects. 

2. For Ay < 5 mag, in the 2 sides models x£ dge rises with Ay. This is physically consistent 
with the fact that cooling photons have a greater and greater path to cross before 
escaping from the far side of the cloud. 
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3. For Ay < 5 mag, in the 1 side models T£ ge is significantly higher at low Ay (50 K, 
about 15%), then drops brutally. This is a numerical effect. One side models make the 
assumption that the cloud extension is large compared to the size of the slab computed. 
However, one still needs to compute some escape probability for cooling lines towards 
the far side. This probability is estimated from level populations in the deepest part 
of the cloud; we use the last few points computed and extrapolate to a large (rather 
arbitrary) depth. If the total depth computed is too small, then the optical depth may 
be overestimated leading to an underestimated cooling. 

Here, C + contributes significantly to the cooling function at the edge, but the size of the C + 
slab is about 3 mag. Photons crossing those 3 mag are essentially free to escape (and thus 
contribute to cooling the gas at the edge), but this (relatively) large escape probability can 
only be computed in models with a total Ay larger than the C + /C transition zone, so that 
the drop in line opacity may be detected. 

As a consequence, the user should remember that it is unwise to compute 1 sided models 
with a too small total Ay, even if the radiation field is large on one side. 

6.2.2. H 2 transitions 

Detailed balance of H 2 is quite complex due to the variety of processes involved. Be- 
sides transitions between ro- vibrational levels, one must include level specific formation and 
destruction processes which have a strong influence on populations. We include: 

• UV pumping to Lyman and Werner bands followed by de-excitation towards either 
the continuum (leading to dissociation) or a ro-vibrational bound state of the ground 
electronic state. 

• Electric quadrupole radiative cascades within the ro-vibrational levels of the ground 
state. 

• Collisional transitions induced by H, H 2 , He, H + , Hjj", including reactive processes 
leading to ortho/para transfer. 

• Level specific formation on grains. 

The H 2 photodissociation rate (used in the chemistry) is a by product of this computation. 
Given the steady state populations of H 2 , it is possible to compute the net energy balance 
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between processes leading to energy loss (typically quadrupole transitions), and energy gain 
(pumping by external UV radiation followed by collisional deexcitation). Thus, radiative 
transitions lead either to heating (mainly close to a brightly illuminated edge) or cooling (in 
deeper parts of the cloud). Note that radiative transitions may arise from H 2 newly formed 
on grain, so that the heating rate from H2 formation above must include both internal and 
translational energy to keep a correct balance. 

6.2.3. Thermal coupling to grains 

We do not compute a detailed energy balance for grains. Dust temperature is evaluated 
using the empirical expression of Burton et al. (1990). Thus, we are not able to compute an 
IR continuum spectrum. However, that temperature is used to compute an energy transfer 
between gas and grains in a manner similar to Burke & Hollenbach (1983). 

7. Chemistry 

We consider mainly elementary two body processes with the exception of possible 3 body 
reactions involving two atomic hydrogen atoms. Reactions are divided into several different 
classes, and read from an input file described in Appendix E. Chemistry for isotopes of C 
and O is automatically generated if 13 C or 18 O are included in the species list. However 
selective reactions involving the main and rare isotopes (leading to isotopic fractionation) 
need to be retrieved and specified. 

From chemical reactions, for each species X, chemical balance leads to: 

at 

with F x the formation rate, and D x the destruction rate. 

With N s species included, built from iV a atoms, Eq (3) leads to a non-linear system of N s 
equations f (n(X)) = to solve. However, the system is under-determined, and only N s — N a 
equations are independent. The iV a missing conditions come from conservation equations of 
each atom. Numerical stability is ensured by excluding at each point the evolution equation 
of the most abundant species for each atom. A Newton-Raphson scheme is used to compute 
the solution. 

This is a dynamical system, whose steady states are solutions of our stationary problem. 
Depending on physical conditions, there may be one or more steady state at each point into 



(3) 
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the cloud as described in Le Bourlot et al. (1995a). Computation starts from the edge of the 
cloud, so that the solution at r = is the "High ionisation Phase" (HIP). As r is increased, 
that solution is followed by continuation even if a new "Low ionisation Phase" (LIP) exists. 
However, deeper into the cloud, there may be a point where the HIP solution ceases to exist 
and only the LIP one is possible. At this critical value of r = r c , the bifurcation from one 
type of solution to the other is characterised by a discontinuity of many variables. Although 
not physical, this is not a bug, nor a weakness of the model. It is a natural consequence of the 
lack of coupling between adjacent regions: it could be solved with the introduction of a small 
diffusive term that would smooth the discontinuity (see Decamp & Le Bourlot (2002) for an 
example of the effect of turbulence on chemistry). However, we feel that the computational 
burden is not worth the trouble. In real clouds both solutions could coexist if, e.g., different 
clumps of gas with different histories were mixed. In that case, reaction-diffusion fronts are 
bound to arise whose dynamics is, to our knowledge, still unknown. 

Once convergence has been reached for radiative transfer, chemical equilibrium and 
thermal balance, the chemical network can be analysed at each point of the cloud using the 
post-processor code. 

The different classes of chemical reactions are described in Appendix E. We focus here 
on two important reactions: formation of H 2 and photodestruction processes. 

7.1. Formation of H 2 

It is customary to write the H 2 formation rate (in cm _3 s _1 ) as _R/n H n(H). The com- 
monly adopted value of Rf, derived from Copernicus observations, is 3 10~ 17 cm 3 s _1 (Jura 
1975). Recent works from FUSE observations (Gry et al. 2002) have confirmed that order of 
magnitude. However, this value results from observational data integrated over whole lines 
of sight. It includes in one single term different physical processes (adsorption of an H onto 
a grain, migration on the surface, reaction with another particle which is not necessarily an 
hydrogen atom, and desorption of the resulting H 2 ). All these processes may vary with depth 
into the cloud as the grain surface certainly does, so the local H 2 formation rate depends on 
r and also on various adopted parameters of the model, particularly the grain model. 

We have chosen not to impose a single constant value of Rf. H 2 production results from 
independent processes: 

H + dust -> H ad (had) 



Had + H ad — > H 2 (k s ) 



-20- 



where "H a( j" means an atom of hydrogen adsorbed on dust, k s is the reaction rate 
(in cm 3 s _1 ) at the surface of the grains and /c a d the adsorption rate (in s -1 ) of H atoms on 
grains. The H 2 abundance may depend also on other processes (various desorption processes, 
or other reactions). If only these two processes are included and steady state applies, then 
the production rate of H2 may be computed by elimination of H adsorbed: 



dt 



= k s n(R ad ) 2 = ^k ad n(H) 

form 



The rate is independent of how adsorbed H atoms eventually manage to reach one 
another. This gives: 

^H 2 1 , 

= -s {an g ) v n n{R) 

form 



dt 



where s is the sticking coefficient of H upon collision with grains, vu its mean velocity, 
and {crn g ) the mean cross section of grains per unit volume. Results from Le Bourlot et al. 
(1995b) show the cross section is given by: 

. . 3 lAm H G 1 

(an g ) = n H 

4 Pg V^min ^ ^max 

where G is the dust to gas mass ratio, p g the density of grain material and the size distribution 
law of grains from Mathis et al. (1977) has been used with a = 3.5. Using G = 10~ 2 , 
p g = 3 gem -3 , a m i n = 3 10~ 7 cm, a max = 3 10 _5 cm one finally obtains: 



dn(R 



2 



dt 



= sl.410 -17 y/T K n H n(H) 

form 



This gives about 410 _17 cm 3 s _1 for Rf at 10 K where s should be close to 1. However, 
it does not include other reaction paths for " H ad " (mainly direct evaporation from warm 
grains) which may lower that rate. Furthermore, at higher temperature the sticking coeffi- 
cient is expected to decrease. The local value of Rf is thus one of the results of the model, 
and depends on the physics and reactions included. 

An effective constant value of Rf may be enforced by using only the two reactions 
above and adjusting the sticking coefficient. For a given grain model, the value of < an g > 
is easily computed and Rf adjusted from the parameter 7 of the adsorption reaction (see 
E.5). However, this procedure should be avoided in normal usage since it is unphysical and 
introduces inconsistencies. 

As an example, Fig. (5) shows the resulting Rf with various sticking coefficients and 
(6) shows the effect of varying s on the H/H 2 transition. The test model has n H = 10 4 cm -3 
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and x — 10. s = 1 leads to a high formation rate, with a H/H 2 transition very close 
to the cloud edge. The "canonical value" of 3 10~ 17 is recovered at the edge by adjusting 
the sticking coefficient to s = 0.213, but this is rather artificial. Our empirical prescription, 
s = \ }? n rr \ gives a constant formation rate close to 3 10~ 17 in the transition region, where 

V max(10,lk) ° ° ' 

most of H 2 emission comes from. Various other laws have been proposed for s, however that 
parameter is still very uncertain (see Hollenbach & Tielens (1999) and references therein). 



7.2. H/H 2 transition 

It is well known that the position of the H/H 2 transition depends on the ratio n^/x 
(Hollenbach & Tielens (1999)). The higher this ratio, the closer to the edge of the cloud is 
the point where / = jv(H)+2iv(H 2 ) = H° wever > in a photon dominated region, the same 
abundance does not translate to the same excitation and in a denser cloud submited to a 
higher radiation field, UV pumping is more efficient for the same column density than in a 
less excited cloud. 

This may be illustrated by computing the H 2 excitation in a small slab of gas of constant 
thickness (say A v = 10~ 2 ) placed at various distances from a star. To be specific, we used 
a late B star ressembling HD 102065 (see also Nehme et al. 2005). runs from 10 3 to 
3 10 6 cm -3 , and d from 5 10" 3 to 1 pc, which gives an equivalent UV field in Draine's unit 
varying from ~ 7000 to 1. Fig. (7) shows the abrupt transition from atomic to molecular 
hydrogen. The fu 2 = 0.5 isocontour shows that the transition is well defined, although not 
a straight line. Fig. (8) shows that to reach the same amount of H 2 on the line of sight, 
density must rise faster than the radiation field. 

Increasing the radiation field also increases radiative pumping of excited rotational lev- 
els of H 2 through cascades. It is thus possible to reach higher excitation temperature for 
a given molecular fraction. Fig. (9) shows iV(H 2 (J = 3)) rising first then decreasing as 
excitation grows. This intermediate level is first populated from lower lying levels, then gets 
depopulated in favour of higher lying ones. This is clearly illustrated in Figure (10) which 
shows the ratio of column densities (J = 5)/(J = 3). Nehme et al. (2005) use these results to 
try and understand observations of rotationally excited H 2 towards HD 102065. They show 
that the observational constraints are met by their models. However, it is unlikely that the 
required high pressure in a slab of gas close to the star will last long enough. Furthermore, 
the observed ortho to para ratio requires that, upon formation on grains, H 2 be released with 
a ratio of 1 instead of 3, which is not supported on physical grounds. 
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7.3. Photodestruction processes 

7.3.1. Dissociation in lines 

For H 2 , CO and their isotopes, photo dissociation proceeds in lines by absorption of 
an ultraviolet photon from the ground state to an electronicaly excited state, followed by 
fluorescence either to vibrationally excited bound states or non radiative transitions in a 
dissociative state via non-adiabatic couplings as described in van Dishoeck (1988). These 
processes are described in details in Abgrall et al. (1992) for H 2 , Lee et al. (1996) for CO 
and its isotopes and Le Petit et al. (2002) for HD. The numbers of discrete levels included 
in a computation for each molecule are parameters of the model. 



7.3.2. Continuous processes 

If ionisation or dissociation cross sections are known, the corresponding photodestruc- 
tion rate is obtained from direct integration over the radiation field by: 

1 f Xt 
k = — / a\ u\ A dX 

h J912 

where k is the destruction rate (s _1 ), a\ the photo-destruction cross section (cm 2 ), u\ the 
radiative energy density (ergcm~ 3 A -1 ) and At the destruction threshold (A), which should 
be larger than the Lyman cut-off. This is done in our code for CI and SI with cross sections 
taken from TopBase (http://vizier.u-strasbg.fr/topbase/xsections.html). The CI ionisation 
cross section happens to be constant (o"c = 1.610~ 17 cm 2 ) from the Lyman limit to the 
ionisation threshold at A t = 1101 . For SI, the section is shown on Fig. (11) (A t = 1168.2 ). 
The resulting ionisation rate for a test model is shown on Fig. (12). 

If the full cross section is not included into the code, then we use the results of van 
Dishoeck (1988). Fig. (12) shows a comparison of the two computations for S for a Ay = 1 
slab of gas illuminated by a standard radiation field. It is seen that the effect is small here. 
However, the full computation includes the effect of protection by H 2 in regions of a cloud 
where that molecule is already abundant and radiation still penetrates between the saturated 
lines. Then it may lead to significant differences in some line intensities as seen on Fig. 3. 
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7.4. Uncertainties in the chemical rates 

Chemical networks are presently available on different websites. The latest UMIST 
database for Astrochemistry is available at http://www.rate99.co.uk/ and described in Le 
Teuff et al. (2000). The chemical reaction rate coefficients are displayed as analytical 
formulae of the temperature as k (J^)" exp (— /3T K ) and are supposed to be relevant in 
a specific temperature range. Chemical networks are also downloadable from the web 
site of prof. E. Herbst from Ohio state University (OSU) at http://www.physics.ohio- 
state.edu/ eric/research. html, with different versions and some comments included. However, 
these data have been initially derived for dark cloud conditions for which very few mea- 
surements are available. Some additions have been implemented regarding neutral-neutral 
reactions involving activation barriers and/or endothermicity (mainly in the UMIST data 
base) and photodissociation probabilities, which are of special interest for PDR models. The 
corresponding photodissociation rates are expressed as k exp (— (3Ay) where k is in s -1 and 
the exponential factor depends on the grain properties, van Dishoeck (1988) describes in 
detail the derivation of the photodissociation and photoionization probabilities and derives 
expressions for different species. Such expressions are very convenient but may be very crude 
as in the case of H2 and CO where the formulae displayed in OSU chemistry are not ad- 
equate. Note that UMIST chemistry does not contain any photodissociation rate for H 2 . 
In our previous studies and the present paper, we compute explicitly the photodissociation 
probabilities of H 2 , HD, CO and the photoionisation probabilities of atomic carbon and 
sulfur as described in Sect. 7.3. 

The chemical network used in the present paper is also downloadable from our website 
and maintained by us. However it is clear that the corresponding data may suffer from many 
uncertainties due mainly to the large temperature range involved in PDRs. The sensitivity of 
the model results to such uncertainties is seldom considered. Roueff et al. (1996) have studied 
the fluctuations in steady state computed abundances resulting from the uncertainties in the 
chemical reaction rate coefficients in the case of dark cold cloud conditions. In such cases, 
the variation of a reaction rate coefficient by a factor less than 2 may even change the 
chemical phase describing steady state, leading to huge differences in the model results. A 
specific example is the dissociative recombination rate coefficient of the Hjj~ molecular ion 
which has been the subject of many theoretical and experimental studies (see for example 
Le Petit & Roueff 2003). In PDRs conditions, we have shown in the case of the horsehead 
PDR (Teyssier et al. 2004), that the results obtained by using the UMIST and OSU 2003 
chemistries are consistent with each other when considering carbon chain abundances. The 
crucial step in PDR is to accurately describe the H, H 2 , C + , C, CO variations which require 
specific treatment as discussed in the present paper. 
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8. Diffuse clouds 



Diffuse and translucent clouds are permeated by the UV interstellar radiation field, 
and are defined as moderately dense media — 10 — 100 cm -3 ) with Ay of a few 1-3. 
Visible and UV observations in absorption, have probed their molecular content: H 2 , HD, 
CO (Copernicus, IUE, FUSE, HST), CH, CH+, CN, C 2 , C 3 , OH, H 2 0, H+ etc. Compared 
to dark clouds, the chemistry of these regions is relatively simple and so they are good places 
to test the physics of the numerical models. 

In order to discuss the trends which can be eventually compared to observations, we 
present a grid of 56 models corresponding to a total visual extinction of 1.0 with different 
densities, nn, and radiation scaling factors, x- These two parameters take respectively the 
values : 20, 50, 100, 200, 500, 1000, 1500 and 2000 cm~ 3 and x =0.1, 0.2, 0.5, 1, 2, 5 and 
10. Thus the ratio nn/x ranges from 2 to 20000. In these low visual extinction conditions, 
the incident radiation field has to be taken impinging on both sides of the cloud. The FGK 
approximation is used (cf Sect. 4.4). This allows us to derive the general trends and to save 
computing time. The elemental abundances and the dust properties are fixed and given in 



T 01 , which is often used as a measure of the kinetic temperature (cf. Rachford et al. 
2002), is defined as: 



We check this assumption by plotting on Fig. 13 the calculated T i as a function of the 
mean kinetic temperature of the gas, T™ 3,11 , derived from thermal balance: 



T i is a good approximation to X£? can for T™ an below 100 K. In the other cases, T i is 
a lower limit as the J = 1 population is not thermalized. They correspond to low ratios 
nn/x an d processes such as radiative pumping in Lyman and Werner bands followed by 
cascades compete with pure collisional excitation. The mean temperature, T 01 , observed by 
Copernicus is 55 ± 8 K and observed by FUSE is 68 ± 15 K (Rachford et al. 2002) and thus, 
assuming T i is close to the kinetic temperature is justified. 



Table 6. 



8.1. T i vs T K 
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8.2. H 2 excitation 



Fig. 14 compares the excitation diagram towards HD 96675 (Gry et al. 2002) with some 
diagrams obtained with our grid of models with Ay = 1 for \ between 0.2 and 2.0. HD 96675 
has been chosen because on this line of sight the probed gas is unperturbed by any special 
radiation field. Among all our models, we select those with densities giving the closest values 
to the observed J = and J = 1 levels. The purpose of this figure is to show that, in classical 
diffuse clouds, radiative pumping is unable to explain the observed excitation of H 2 in its 
levels J > 3. Other mechanical processes such as C-shock (Flower & Pineau des Forets 1998) 
or turbulence (Joulain et al. 1998; Casu 2003) are required to explain it. 



The atomic and molecular abundances of hydrogen are the result of the balance between 
the formation of H 2 on dust and photodissociation. As a consequence we expect that their 
column densities are a function of the ratio of the density to the scaling factor of the incident 
radiation field: n H /x- Observers define the molecular fraction as : 



Fig. 15 displays the variation of the molecular fraction as a function of n-^/x- First, as 
expected, / reaches a value of 1 for high nu/x ratio. From FUSE and Copernicus observations 
of diffuse and translucent clouds, Rachford et al. (2002) brought to the fore a group a 10 
lines of sight 2 with / ~ 0.7. According to Fig. 15 this corresponds to n-^/x — 60. So if 
we assume that x — 1? this gives riu around 60 cm -3 which is the order of magnitude of the 
expected value of the density in diffuse/translucent clouds. Secondly, we see that models 
with different parameters but the same ratio nn/x gi ye the same results: for a given visual 
extinction, the molecular fraction of a cloud is controlled by its ratio n H /x- 



Fig. 16, 17, 18, 19 and 20 present the column densities of some observed species in 
diffuse clouds as a function of nu/x- The column densities of CH, CN, C and C + are 



2 These lines of sight are C Oph, o Per, C Per, HD 24534, HD 27778, HD 62542, HD 73882, HD 96675, 
HD 154368 and HD 210121. 



8.3. Atomic to molecular fraction 




8.4. Chemical results 
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directly controlled by the ratio riu/x- This behaviour appears when species are formed by 
two body reactions and are destroyed by photo-processes (or the opposite as for C + ). For 
species such as CO, OH and NH there is no simple relationship since other mechanisms, 
which do not involve two body reactions but specific reactions such as cosmic ray ionisation 
(for OH and NH), come into play . The case of CO is particular and can be explained by 
its forming reaction. For models with nu/x — 100, its formation occurs via OH + C + . 
Whereas for higher values (n H /x — 10 4 ) the reaction O + C 2 is the main route of formation 
of CO. Thus the behaviour of iV(CO) vs n H /x follows the one of OH or C 2 depending on 
the value of nu/x- 

8.5. UV radiative transfer in H 2 lines 

We checked the validity of the FGK approximation for the model nn = 100 cm~ 3 and 
X — 1, Ay — 1, by solving the full radiative transfer as described in Sect. 4 up to J = 5. The 
photodissociation probabilities of H 2 , HD, CO and photoionisation probabilities of C and S 
are the same at the edge of the cloud (Tab. 7) in the two treatments. However, the values 
differ significantly at Ay greater than 0.05 as displayed on Figs. 21 and 22. 

Table 8 presents the column densities obtained in the two cases. The column density of 
CO is increased by a factor of 2 when proper account of the radiative transfer is performed. 
This produces also a decrease of iV(H) and an increase of N(R 2 ) and N(C). For molecules 
(CH, OH, etc.) whose photodissociation rate is given by the formula ■ye~ f3Av ((van Dishoeck 
1988)) column densities are obviously similar with the two treatments. 

Extending our photodestruction cross-sections database to compute photorates by direct 
integration over the true radiation field is under way. 

8.6. One side and Two sides models 

The numerical model allows the incident radiation field to come from one side or from 
two sides of the plane parallel cloud. In this latter case, the computing time is increased 
significantly in order to converge. A two sides model is required when the radiation field 
coming from the other side is still significant, i.e. for low value of Ay. We study this effect 
by comparing two models (tir = 100 cm -3 , and x = 1) : model 1 where the radiation field 
comes from one side and model 2 with the radiation field coming from both sides of a cloud 
with a visual extinction Ay*. In order to compare column densities between the two models, 
in model 1, they are computed integrating abundances up to Ay f /2 and are then multiplied 
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by2. 

Results (Table 9) show that differences are important for iV(H 2 ) only for low visual 
extinction clouds (< 0.5). However, for other molecules such as CO, the differences are 
significant at much higher Ay (we have a factor of 2 differencies in column densities at 
Ay = 1). For diffuse lines of sight, using a 2 sided model is mandatory. 

9. Dense PDRs 

In this section, we consider dense clouds in a young stellar environment. We present a 
grid of models with radiation field coming from one side to model a dark cloud with high 
extinction illuminated by a young star. The total visual extinction is chosen as 20. Adopted 
densities are n H = 10 3 , 3 x 10 3 , 7 x 10 3 , 10 4 , 3 x 10 4 , 7 x 10 4 , 10 5 , 3 x 10 5 , 7 x 10 5 , 10 6 , 
3 x 10 6 , 7 x 10 6 and 10 7 cm -3 . The adopted multiplicative factors of the Draine radiation 
field are x = 10 3 , 3 x 10 3 , 7 x 10 3 , 10 4 , 3 x 10 4 , 7 x 10 4 , 10 5 , 3 x 10 5 , 7 x 10 5 , 10 6 , 3 x 10 6 , 
7 x 10 6 and 10 . These parameters correspond to a wide range of physical conditions. Other 
parameters are given in Table 6. 

9.1. Temperature at the edge of the clouds 

At the edge of the cloud, the main heating process is the photoelectric effect on dust. 
The kinetic temperature at the edge of the cloud obtained in the grid of models is presented in 
Fig. 23 and varies between 100 and 2500 K. It is often used as an indicator of the excitation 
of the medium. However, as shown in Fig. 24, the temperature maximum is not located at 
the edge of a cloud but rather a bit deepper. For highly illuminated clouds, grain ionization 
can be so large at the edge that the photoelectric effect is less efficient than deeper where 
the charge of the grains decreases as displayed in Fig. 25. 

9.2. Line intensities 

We compute the intensities (in ergs cm -2 s" 1 sr _1 ) of the cooling transitions assuming 
that the PDR is seen "face on". The results for C + , C, O, CO are in good agreement with 
the results of Kaufman et al. (1999) and are not reported here. We focus on the molecular 
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hydrogen infrared transitions over a wide range of densities and UV enhancement factors 3 . 

The 1-0 S(l) transition at 2.12 /im has been widely studied in a variety of environments 
and the corresponding intensity is displayed in Fig. 26. 

It is often said that the ratio of the intensities of the H2 lines 2-1 S(l) on 1-0 S(l) is a 
good candidate to discriminate between PDR and shocks: a value of 0.11 is assumed to be 
characteristic of shocks (Kwan 1977) and a value of 0.54 is appropriate for radiative pumping 
(Black & van Dishoeck 1987). Fig. 27 presents this ratio in the plane n H — X obtained with 
the grid. For n H < 10 5 the ratio is almost constant at a value close to 0.56. The ratio is 
decreasing at higher densities with increasing radiation fields and a ratio of 0.1 is obtained 
for rin > 10 5 cm" 3 and x > 10 5 . So, if such parameters cannot be justified, other physical 
processes such as shocks are required to explain such a ratio. 

Figs. 28, 29, 30, 31 and 32 present respectively the ratio of the intensities of the 1-0 
S(0), 1-0 S(2), 1-0 S(3), 2-1 S(2) and 2-1 S(3) H 2 lines to the 1-0 S(l) H 2 line which may 
be detected in good observing sites. These ratios are almost always smaller than 1. The 
sensitivity of the ratios to density is moderate below 10 5 cm" 3 but more significant at higher 
densities. The 1-0 S(2) / 1-0 S(l) and 1-0 S(3) / 1-0 S(l) ratios could be used for high 
density determinations. On the other hand, these ratios are highly dependent on x an d the 
observation of several of these transitions can lead to a good guess of the intensity of the 
incident radiation field. Ratios involving para and ortho transitions are different for each 
case. 



9.3. Ratios of antenna temperatures 

CO is often used as a tracer of H 2 and deserves special discussion. We display the 
isocontours of the intensity of the CO (2-1) transition in Fig. 33 in the n H - X plane. We 
checked that the results are in agreement with Kaufman et al. (1999). We also display the 
ratio of the antenna temperatures of 3-2 and 6-5 transitions to the 2-1 transition in Fig. 34 
and 35. 

The ratios of the antenna temperatures (Ta) are defined as (for example in the case of 
Tj^ 1 /T^~°) , using the Rayleigh- Jeans approximation: 



3 The referee mentioned a forthcoming paper by Kaufman et al. (2006) in which mainly rotationnal 
transitions of H2 are displayed. In the one case in which there is an overlap over a limited range of parameters, 
the 1-0 S(l) line intensity, the agreement is quite good. 
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T\-° U-J h-o 

where vi-j is the frequency of the i — > j transition and Ii-j is the intensity of this same 
transition. 

In the hypothesis of an homogeneous and optically thin medium, the ratio of the pop- 
ulations is given by : 

n2 ^ T^ 1 Ai_ i/ 2 -i 

ri! ~ A 2 _x I/!_ 

Whereas the T^~ 2 to T^ 1 ratio is only slightly dependent on the density and illumination 
conditions, much larger variations are obtained for the highly excited (6-5) to (2-1) antenna 
temperature ratios. Thus highly excited CO transitions are valuable tests of dense PDR 
conditions. 



10. Conclusion 

We have presented and discussed the new features implemented in the Meudon PDR 
code first described in Le Bourlot et al. (1993). The principal highlights are: 

• Treatment of the UV radiative transfer : a 2 sides illumination of the gas slab is in- 
cluded as well as the possibility of an additional UV source impinging perpendicularly 
to the surface of the "cloud". The treatment of the UV radiative transfer includes the 
effect of discrete transitions of H and H2 in the computation of the radiative energy 
density. This possibility, which is computing time intensive, is optional and is tuned by 
the number of rotational levels of H2 explicitly involved. The approximate treatment 
following the FGK approximation allows for rapid computation with a reasonable ap- 
proximation to the main features of the UV radiative transfer involving self-shielding 
effects but neglecting discrete line overlaps between H, H 2 and CO. Photodissociation 
of H2, HD, CO and its isotopes by summation over all relevant transitions, as well as 
photoionization of C and S from the integration of photoionization cross-sections over 
the UV radiative field intensity are then obtained. 

• Grain properties : Scattering properties of grains are explicitly introduced for carbon 
and silicate particles for sizes ranging from the nanometer to the millimeter region. 
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Corresponding data are taken from Weingartner & Draine (2001). PAH properties are 
mimicked by the nanometric carbon aggregates data. The photoelectric effect as well 
as the grain charge distribution are then computed from the actual value of the UV 
radiation field without analytic formula. 

• Chemical processes : Chemical processes involving grain surfaces are included as a 
function of the grain size distribution. In this respect, H 2 formation results from 
subsequent adsorption of H on the grain surface and reaction between the two adsorbed 
atoms. Three body collisions are implemented for dense regions where two hydrogen 
atoms and a third body may contribute to chemical formation . 

• Thermal balance : An effort to homogenize the various cooling processes has been 
achieved and includes the excitation and subsequent emission of forbidden visible tran- 
sitions of the main atomic and ionic constituents. 



Examples relevant both for diffuse and dense cloud conditions are presented. In the 
diffuse cloud conditions, we discuss possible chemical diagnostics as a function of nn/x- We 
show that whereas the molecular fraction, fractional column densities of C + , C, CH, CN 
scale smoothly with n H / x, no specific trend is found for CO and OH. We have displayed 
emissivity ratios of rovibrational transitions of H 2 in a parameter space n#, x, relevant to 
dense PDRs. These plots are a guidance for the interpretation of observations. However, 
we recommend our readers to rather run their own model with the proper parameters of the 
studied line of sight. Our code is indeed downloadable at http://aristote.obspm.fr/MIS/ 

Several improvements are still in progress. These concern the treatment of photodisso- 
ciation processes with the inclusion of appropriate cross sections when available. This allows 
us to include the effect of grain extinction properties without any approximation. The role 
of infrared pumping is also under study together with a consistent treatment of the grain 
properties and gas phase elemental abundances. Finally, we also plan to consider the depen- 
dence of the grain size distribution on the space position within the cloud: there is increasing 
observational evidence that very small grains and PAH are abundant in the illuminated parts 
of the cloud whereas big grains are more likely present in the shielded regions. We neverthe- 
less consider that the present PDR code is a useful tool to the scientific preparation of future 
spatial missions such as Herschel and Spitzer and look forward to incorporate it within the 
Virtual Observatories. 
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A. 



Cloud structure 



The cloud is ID, plane parallel, with sharp edges. It may be semi-infinite or have a 
finite extent Ay 1 in magnitude, in which case radiation comes from both sides. Geometry 
and sign convention are displayed on Fig. (1). The origin is at the left edge of the cloud. The 
observer is always on the negative side. Integrated line intensities are computed following 
an angle 9 from the normal of the surface of the cloud. We define x + an d X~ respectively 
as the enhancement factors of the Draine radiation field on the positive and negative sides 
of the cloud. If x + = 0, the cloud is semi-infinite and Ay 1 is the limit up to which the 
structure is computed. If x + 7^ 0, Ay 1 is the cloud total visible extinction. It is also possible 
to introduce an illuminating star at a distance from the cloud. The star is described by 
its spectral type, corresponding to an average star radius and blackbody temperature. It 
may be either on the same side as the observer (d* < 0) or the opposite side (d* > 0). No 
star is set by d* = 0. 

The gas column density follows from the relations Cd = and Ry = , where 

A^h is the total hydrogen column density in cm -2 (iVn = -W(H) + 2 N(H 2 )). The standard 
galactic values are Cd = 5.8 10 21 cm~ 2 mag -1 (Bohlin et al. 1978; Rachford et al. 2002) and 
Ry = 3.1. The relation between the optical depth and the path length is: 



where nn(rv), in cm 3 , is the total hydrogen density at a visible optical depth r v . 



The steady state solution is reached only after a number of iterations over the whole 
structure of the cloud. This comes from the necessity to know in advance the optical depth 
in lines towards both sides of the cloud in order to compute the thermal balance. This in 
turn requires knowledge of column densities computed from abundances. Hence some kind 
of initialisation is required, followed by iterations. Initialisation is rather arbitrary. 

Convergence of such a process is not guaranteed. However, one finds that most physical 
quantities do converge in very few iterations towards a value close to the final one so that 
most numerical problems arise during the first or second iteration and are often cured by 
skillful tuning of the initial guess. 

One quantity however converges significantly slower than any other and thus controls 
the number of iterations: it is the position of the H/H 2 transition which is sensitive to state 




B. Numerical convergence 
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dependent photo-dissociation of H 2 . Self shielding in lines depends on J with a protection 
significantly more efficient for ortho-H 2 compared to para-H 2 . The consequence is a peak in 
the ortho to para ratio occurring just after the transition from atomic to molecular hydrogen, 
as already pointed out by Abgrall et al. (1992). Iteration after iteration, that transition is 
pushed further into the cloud as level column densities are computed more accurately. The 
situation is illustrated on Fig. 36 for a cloud model with % = 10 4 cm~ 3 and x — 10 4 . 20 
iterations were required to reach a satisfactory stability. Less stringent physical conditions 
(i.e. usually a lower radiation field) lead to convergence in less than 10 iterations. 



C. FUV radiation field 

The basic physical quantity for radiation is the specific intensity /(A), measured in 
erg cm -2 s" 1 sr" 1 A -1 . Throughout most of the code, we use wavelengths expressed in A as 
the independent variable. The most useful derived quantities are the mean intensity J(A) 
(in erg cm -2 s" 1 A -1 ) defined by 

J(A) = i- J J(A)dfi 

and the energy density u(X) (in erg cm" 3 A -1 ) defined by 

u(\) = - f I(X)dn = —J(X) 
c J c 

If (and only if) the radiation field is isotropic, I = J. /(A) is the quantity to use to 
solve the radiative transfer equation. However, most physical quantities depending on the 
radiation field can be computed from u(X). 

Our standard impinging radiation field is the value given by Draine (1978). Various 
incompatible expressions are found for that field, as discussed in Kopp (1996). We have 
chosen to use the formula of Sternberg & Dalgarno (1995) which reads 

1 /6.30010 7 1.0237 10 11 4.0812 10 13 
J (A) = - — ^ r , + 



4tt V A 4 A 5 A 6 

where A is in A, and J(A) in erg cm" 2 s" 1 sr" 1 A -1 . This expression is used from the Lyman 
cut-off up to the limit given by Draine (1978) at 2000 A. Longward, we use 



/(A) = 1.38243 10" 5 A" 3 
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This field may be scaled by a factor x, an d an additional radiation from a nearby star 
taken as a diluted Black Body can also be introduced. The star is characterised by its 
effective temperature T ef[ , its radius and its distance to the cloud surface d*. 

It is common to measure the strength of the radiation field in reference to Habing (1968). 
We take the Habing standard value as 5.6 10~ 14 ergcm~ 3 between 912 A and 2400 A. Note 
that this cut-off corresponds to 5.166 eV, whereas it is often given as 6eV. Although small, 
this difference is one amongst the many inconsistencies between various codes. Here, we 
define the radiation scaling factor at some point by 

i /-2400 

G = - / u(\)d\ 

5.6 10- 14 J 912 1 ; 

The "usual" Go parameter is taken as the value of G without cloud (in free space). Since u(\) 
comes from an integration over 47rsr, this value differs from G at the cloud surface where 
half of the radiation is partly screened by the cloud itself. For a semi-infinite cloud and an 
isotropic impinging radiation field, taking into account back scattering of the radiation by 
dust, G sur f acc is usually close to 0.54 Go (depending on the dust characteristics). This may 
lead to (close to) a factor of 2 of difference in the energy input of otherwise seemingly similar 
models, and is a major source of discrepancies. 



D. Spherical Harmonics solution of the transfer equation 

We follow Roberge (1983) for most of the development, but take into account the fact 
that grain and/or gas extinction and scattering properties may vary with position into the 
cloud. Here, we do not take into account the possibility of embedded sources. The following 
development applies to a plane parallel cloud with coherent scattering so that the dependence 
on wavelength is omitted. In the radiative transfer equations, Eq (1) and Eq (2), u){r) = 
K o + a K D +a D is the troublesome contribution. This effective albedo includes absorption by the 
gas in lines, and is thus smaller than the usual dust albedo uj d {t) = K D+ a D ■ Boundary 
conditions are (with r max the total slab optical depth at A): 

f I(t = 0,//) = fi<0 
If Tmax = oo, I + = 0. We expand I and p(/x, //) in Legendre polynomials Pj by: 

oo 

J(r, Ai ) = ^(2/ + l)/ i (r)P(/i) 
1=0 
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1=0 

Flannery et al. (1980); Roberge (1983) show that this leads to a infinite set of equations 
whose general term is: 

I fU + (/ + = (21 + 1)(1 - W (r)(7,) /, 

Our expression differs from Roberge (1983) (Eq 9) by the explicit r dependence of u, so that 
this is not a constant coefficient equation. We now stop the expansion at an odd value L, 
and write it as a linear system: 

A(t) f = f 

From there, Roberge (1983) (Eq 11, 13, 14, 15, 16) hold, with all coefficients hi, k m and Ri m 
now depending explicitly on r. We still have k_ m = —k m and Ri- m = (—l) l Ri tm . Letting 
f = Ry, we have f = Ry' + R' y. The derivative of R introduces a new term in the equation, 
so that Roberge (1983) (Eq 16) is now: 

y' = k(r) y + Q 

where Q = R _1 R'y does not decouple. Each component q m (j) is a known linear combination 
of Umij). This is a linear equation, whose homogeneous part has a solution: 

yjj) = C m exp (^j k m (t) dtj 

where the 2M constants r m may still be chosen arbitrarily and C m = y m (j m ). 
Adding the particular solution, one has: 



y m (r) = exp {j^ k m (t) dtj x 

[C m + f Tm exp (- j s Tm k rn (t) dtj q m (s) ds j 



The integral involving q m is large only if both R' and q m are large, which occurs only 
where lines dominates the absorption and are not yet saturated. A first order approximation 
is obtained by considering this integral as a constant. Iterative refinements are possible 
afterwards, but are not taken into account here. A more accurate treatment will be the 
subject of another paper. 

One sees that, provided the albedo oj(t) is known, a complete solution is still com- 
putable, but requires numerical integrations. Those expressions may thus be used in an 
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iterative scheme. Back to the original variables, we have now: 

+M 



fi( T ) = Yl Rim( r )C m exp( k m (t)dtj 



rn=-M V% 

m = being omitted from the sum. 

To avoid numerical instabilities, we also write: 

/iW=/, + (r) + (-l) , /r(r) 

with 



/i + 00 = ZZ=i Rim{r) C rn exp (- k m (t) dt) 
fr( T ) = 2=i Rim{r) C. m exp (- £ k m (t) dt) 



where we have chosen r m = if fi m < 0, and r m = r max if fi m > 0. 
Constants C m are determined from boundary conditions by: 

M 



^ C m B im — Qi (Dl) 



m=-M 



with 5 im given by: 

• Hi < ; m < 0: 

• /ij < ; m > 0: 



J2(2l + 1) fl(^) R lm (0) 



L / fTmax 

V(2Z + 1) ^ m (0) exp - / fc^df 



/ij > ; m < 0: 



L / /"T max 

J^(2/ + 1) P t (fii) i? im (r max ) exp ^- y 



/ii > ; m > 0: 

L 

^niax / 
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and 



/ (0,/ii) fii < 



In the following, only the mean intensity Ja is needed. This is given by /o which depends 
on i?o,m = 1- So all coefficients Ri m {r) need not be kept during computation. One only needs 



Species (atoms and molecules), initial abundances and chemical reactions to use in a run 
are given in a single file. The first part gives the list of the species with their name, atomic 
composition, initial abundance and enthalpy of formation in kcal mol" 1 . These enthalpies are 
used in the thermal balance to compute the variations of enthalpy of the chemical reactions. 

The second part of the chemistry file is a list of the chemical reactions. Chemical rates 
are computed from the three parameters 7, a, (5 given for each reaction. Photo processes 
computed by direct integration over the radiation field (i.e. photo-dissociation of H2, HD, 
CO and its isotopes and photo-ionisation of C and S) must not be mentioned in this file. 
Reactions on grains may be included, following Le Bourlot et al. (1995b). Most of the listed 
reactions follow an Arrhenius law : 



the coefficients C m and the integrated eigenvalues JJ" k m (t) dt. 



E. Chemical processes 




with k the chemical rate at a point of the cloud at the temperature Tk. 



Here, we would just like to mention how we deal with some specific reactions. 



E.l. Secondary photon processes 



X + hz/ 



products 




-1 



These photons are created deep into the cloud by electronic cascades of H 2 following 
excitation by electrons produced by cosmic rays as described first by Prasad & Tarafdar 
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(1983). Temperature dependence is needed for CO only. Rates are taken from Gredel et al. 
(1989) where an albedo of 0.5 is assumed. 



E.2. Radiative association 



X + Y -> XY + hv 



fc = 7 exp(-/3/r K )cm 3 s- 1 

These reactions are singled out from ordinary gas phase because it is not possible (due 
to the escaping photons) to compute their contribution to thermal balance from simple 
thermodynamic considerations. 



E.3. Endothermal reactions with H 2 

X + H 2 — * products 

(rp \ Oi I max 
^^expH/^-^/r^cm^- 1 

where n\ is the relative population of level / of H 2 , E\ its energy, and / max is the highest level 
of H 2 such that f3 — Ei > 0. Those reactions make the implicit hypotheses that all internal 
energy of H 2 may be used to overcome an activation barrier or an endothermicity. 

E.4. Photoreactions 

X + hp — > products 
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These reactions are used only for species whose destruction rate is not computed directly 
by integration over the local radiation field, x^ are the scaling factors of the radiation field 
with respect to that of Draine on the left and right side of the cloud respectively (see 
Appendix C and A for notations), and Ay ax = 2.5 log 10 er max is the total cloud extinction. 
For a semi-infinite cloud, x + = 0- Note that no factor of \ is needed to take into account 
the fact that photons come from only a half space at each edge, since this is already taken 
into account in the computation of x ± - 



E.5. Adsorption on grains 



X + dust -> X : 



k — 7 (<Jn g ) v 

where X: is a species adsorbed on grains, < an g > is the mean grain cross section per unit 
volume, and v the mean particle velocity. Microphysics (including the sticking coefficient) 
is "hidden" in the term 7, which is read from the chemical input file for most species. We 
use specific expressions for H and H 2 : 



• For H, the default sticking coefficient is 7 = y sup (io x K ) • As f oc V^k~, this leads to a 
constant rate at high temperature. The consequences are discussed in Sect. 7.1. 

• For H 2 , the rate is tuned to give at most a single mono-layer on the grain, as described 
in Le Bourlot (2000). 



E.6. Grain surface reactions 



X : + Y : — > products 



Gho P (X £ho P (Y:)) 

where F r is the fraction of sites on the outermost layer occupied by Y, X is either H or D, 
and t hop is the "hopping" time of the particle from one site to the next. See Le Bourlot et 
al. (1995b, Appendix C) for details and the expression of F r . We use tho P (H) = 2 10~ n s. 



Note that this expression supposes that the time to reach a specific position varies as 
the distance to it and not as the square, as would occur for a random walk on an infinite 
surface. It is also not valid for reactions other than hydrogenation. 



E.7. Thermal evaporation from grains 



X : 



X 



k — 7 




) 




where (3 is the adsorption binding energy of X on the surface, y ^ a vibration frequency 
and T d the dust temperature, see Cazaux & Tielens (2003) for details. 

We thank the referee for his careful reading which improved the english considerably. 
The authors thank Pierre Hily-Blant for numerous discussions and heavy testing of the code. 
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Table 1. Model variables defined or calculated at each point in the cloud. 



variable 


Unit 


Name 


Comment 


T K 


K 


Kinetic temperature 


Variable or parameter 


n H 


cm -3 


Density (Hydrogen nuclei) 


Variable or parameter 


n(X) 


cm -3 


Abundance of X 




ni(X) 


none 


Population of level i of X 


E^(x) = i 


«(A) 


ergcm~ 3 A -1 


Radiative energy density 




a & r 


erg cm -3 s _1 


Heating and cooling rates 


See Sect. 6 
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Table 2. Astrophysical parameters : model definition. 



Parameter Units Name or definition Comment 



X 'Draine' FUV radiation strength see Appendix C 

£ 10~ 17 s _1 Cosmic ray ionisation rate 

Ay mag Extinction Total cloud depth 

fturb cms -1 Turbulent Velocity Should be "local" 

5x none Depletion of atom X Relative to H 

P K cm -3 Thermal pressure Defined by P = n T K 

with n = n(H) + n(H 2 ) + n(He) 



Table 3. Atomic and molecular parameters. 



Constant 


Units 


Name 


Comment 


kxy(T K ) 


s —1 
cm s 


Chemical reaction rate coefficient 


Elementary two body process 


<j x (E) 


2 

cm 


Cross section 


Photoionisation, dissociation, etc... 


Aij 


s- 1 


spontaneous emission transition probability 


Derived: B ih B jU fa 


Qij(T K ) 


cm 3 s _1 


Collisional excitation rate coefficients 


Derived: 
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Table 4. Grain parameters. 



Parameter 


Units 


Name or definition 


Comment 


Typical value 


^min 


cm 


Lower size cut-off 


MRN 


3f0~ 7 


"max 


cm 


Upper size cut-off 


MRN 


3f0~ 5 


d 


cm 


Mean distance between 


Le Bourlot 


2.6 fO" 8 






adsorption sites 






a 


none 


index 


MRN 


-3.5 




none 


dust albedo 


fixed a , M96 


0.42 


9 


none 


anisotropy factor < cos 6 > 


fixed a , WD 


0.6 


G r 


none 


-^grain /-^gas 


fixed a 


O.Of 


Pgr 


gem -3 


grain volumic mass 


fixed 3. 


3 


Ry 




A v /E B -v 


See Appendix A 


3.1 




cm 2 mag -1 


A^h/Eb-v 


See Appendix A 


5.8 fO 21 


C3, 7, yo 


depend 


Extinction curve, 2200 bump 


FM 


Galactic 


<5abs(a, A) 


none 


absorption coefficient 


DL 





References. — MRN: Mathis et al. (1977), FM: Fitzpatrick & Massa (1986, 1988, 1990), Le 
Bourlot: Le Bourlot et al. (1995b), DL: Draine & Lee (1984) upgraded by Laor & Draine (1993), 
see http://www.astro.princeton.edu/~draine/dust/dust/diel.html, M96: Mathis (1996), WD: 
Weingartner & Draine (2001). 

Parameters quoted "fixed" are in fact dependent on grain composition. Modifications to 
the program for a more physical modelization are in progress. 
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Table 5. Collisional processes included. 





1 Dl TO 1 O 

ije veis 


u 
j j 


1 It 


1 1 


n 2 


6 


c 


3 p 3p 3p In 1 Q 
JO, .Ti, 12, J_y2; »J0 


(1) 


(2) 


(3) 


(4) 


(5) (6) 





3p 3p 3p In 1 c 
"li -Tl, -TO, i^2 5 »J0 


(1)(7) 


(8) 


(9) 


(10) 


(9)(H) 


s 


3p 3p 3p In 10 
-<2, -Tl) M)) -^2) »J0 


(1)(+) 


(8)(+) 


(9)(+) 


(io)(+) 


(9)(10)(+) 


Si 


3p 3p 3p In lo 


(l)(t) 


(2)(t) 


(3)(t) 


(4)(t) 


(5)(6)(t) 


c+ 


2 Pl/2, 2 P?,/2i 4 Pl/2, 4 P3/2, 4 Pb/2 


(12) 






(13) 


(14)(15) 


N+ 


3p 3p 3d In 1 c 
■M)) -ri; -T2, -^2, ^0 










(6) 


Si+ 


2 -Pl/2, 2 -P3/2, 4 -Pl/2 5 4 -P3/2 5 4 -P5/2 


(16) 






(13)(t) 


(17) 


HCO+ 


J = -> 20 








(18) 


(19) 


CO 


J = -> 30 


(20) 


(21) 




(22) 




cs 


J = -> 20 








(23) 


(24) 


H 2 


J = -> 29, v = -> 14 a 


(25) 


(25) 




(25) 




HD 


J = -> 9 


(26) 


(26) 




(26) 





Note. — A "-" indicates that no collision is considered. (+): Rates for S are taken from 
O. (f): Rates for Si are taken from C, and some rates for Si + are taken from C + . 

a All ro- vibrational levels of H 2 ground electronic state may be included. Usually they are 
included up to a highest level (l ) chosen on physical ground. 

References. - (1) Launay & Roueff (1977a), (2) Lavendy et al. (1991); Staemmler & 
Flower (1991); Le Picard et al. (2002), (3) Roueff & Le Bourlot (1990), (4) Schroder et al. 
(1991), (5) Pequignot & Aldrovandi (1976), (6) Mendoza (1983), (7) Federman & Shipsey 
(1983), (8) Monteiro & Flower (1988), (9) Chambaud et al. (1980); Pequignot (1990), 
(10) Jaquet et al. (1992), (11) Bell et al. (1998), (12) Launay & Roueff (1977b), (13) Flower 
& Launay (1977), (14) Lennon et al. (1985), (15) Wilson & Bell (2002), (16) Roueff (1990), 
(17) Dufton & Kingston (1991), (18) Flower (1999), (19) Faure & Tennyson (2001); Neufeld 
k Dalgarno (1989), (20) Balakrishnan et al. (2002), (21) Cecchi-Pestellini et al. (2002), 
(22) Flower (2001), (23) Turner et al. (1992), (24) Dickinson et al. (1977), (25) Le Bourlot 
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et al. (1999) and ref. therein, (26) Flower et al. (2000) and ref. therein 
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Table 6. Adopted gas phased abundances relative to % and adopted parameters in the 

grid of models (see also Table 4). 



Parameter 


Value 


He 


0.1 


C(°) 


1.3 (-4) 


N (b) 


7.5 (-5) 




3.2 (-4) 


sw 


1.9 (-5) 


Fe (a) 


1.5 (-8) 


C(s- r ) 


5.0 (-17) 


bikms- 1 ) 


2.0 



Note. — Ref: (a) Sav- 
age & Sembach (1996), 

(b) Meyer et al. (1997), 

(c) Meyer et al. (1998) 
Figures in parentheses 
are powers of ten. 
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Table 7. Edge destruction probabilities in s 



H 2 HD CO C S 



4.2 (-11) 2.6 (-11) 1.1 (-10) 1.7 (-10) 5.1 (-10) 



Note. — Model parameters: n H = 100 cm 3 , x = 1, 
A v = 1. Values in parentheses are powers of 10. 
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Table 8. comparison of FGK approximation and exact radiative transfer. 



FGK Exact FGK Exact 



H 3.5(20) 2.4(20) C+ 2.4(17) 2.4(17) 

H 2 7.6(20) 8.1(20) C 1.0(15) 1.6(15) 

T i 65 62 CO 4.6(13) 9.3(13) 

CH 3.0(12) 3.2(12) OH 2.9(13) 2.5(13) 



Note. - Column densities in cm" 2 and T 01 
in Kelvin obtained with the FGK approxima- 
tion and the exact radiative transfer calculation 
up to J = 5 (see Sect. 4). For both models, 
nn = 100 cm -3 , x = 1, radiation field from 
both sides, Ay = 1 and the parameters in Ta- 
ble 6 are used. 
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Table 9. One and two sides results. 



A tot 






0.2 


0.5 


1.0 


5.0 


7.0 


H 


1 


1 


7(20) 


2.6(20) 


3.2(20) 


4.8(20) 


5.5(20) 




2 


1 


9(20) 


2.9(20) 


3.5(20) 


5.0(20) 


5.7(20) 


H 2 


1 


1 


0(20) 


3.4(20) 


7.7(20) 


4.4(21) 


6.3(21) 




2 


9 


0(19) 


3.2(20) 


7.6(20) 


4.4(21) 


6.3(21) 


C+ 


1 


3 


0(16) 


1.3(17) 


2.4(17) 


8.4(17) 


8.8(17) 




2 


4 


9(16) 


1.2(17) 


2.5(17) 


8.5(17) 


8.8(17) 


c 


1 


9 


2(13) 


3.7(14) 


1.5(15) 


3.7(17) 


8.0(17) 




2 


4 


8(13) 


2.1(14) 


1.0(15) 


3.7(17) 


8.0(17) 


CO 


1 


8 


8(12) 


3.0(13) 


8.6(13) 


2.1(16) 


5.4(16) 




2 


2 


7(12) 


1.3(13) 


4.6(13) 


1.9(16) 


5.2(16) 


CH 


1 


4 


4(11) 


1.8(12) 


4.8(12) 


6.4(13) 


1.2(14) 




2 


1 


5(11) 


8.3(11) 


3.0(12) 


5.6(13) 


9.9(13 



Note. - Comparison of column densities for 1 side 
and 2 sides models. Figures in parentheses correspond to 
powers of 10. For all models, nn = 100 cm -3 and x = 1. 
The second column gives the number of sides of the model. 
For the one side model, column densities correspond to 2 
times the value at Ay/2. 
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Fig. 1. 



Structure and geometry conventions. 
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Fig. 2. — Radiative energy density inside a 0.5 mag clump, riu = 100 cm 3 , irradiated by the 
standard ISRF on both sides. 
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Fig. 6. — Resulting H/H2 transition for various sticking coefficient s. Same parameters as 
Figure 5. Note the large shift in the H/H2 transition. 
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Fig. 7. — H/H2 transition for a A v = 10 2 slab of gas as a function of density rin and distance 
d (pc) to the star HD 102065. The solid line in the plane ci-% corresponds to /h 2 = 0.5. 
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Fig. 8. — nu/x ratio for the iso-/H 2 — 0-5 curve. Fig. 7 shows that for each value of n H 
(lower x axis) a single radiation field leads to /h 2 = 0.5. The associated distance to the star 
is shown on the upper x axis. 
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Fig. 9. — Column density of H2 (J = 3) as a function of for various distances to 
HD 102065. 
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Fig. 13. — Excitation temperature of levels J = 0, 1 of H 2 versus the mean kinetic tempera- 
ture. Each point corresponds to a model (Ay = 1). Densities are 100, 200, 500, 1000, 1500 
and 2000 cm" 3 . The straight line displays T 01 = T^ oan . 
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Fig. 14. — Comparison of the excitation diagram towards HD 96675 (Gry et al. 2002) with 
some models from our grid. 
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Fig. 15. — Molecular fraction as a function of n^/x- Each point corresponds to a model 
(Ay = 1). Densities corresponding to these models are 20, 50, 100, 200, 500, 1000, 1500 and 
2000 cm- 3 . 
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Fig. 18. — Column densities of CN as a function of rin/x- 
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Fig. 19. — Column densities of CO as a function of n^/x- 
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Fig. 20. — Column densities of OH as a function of tih/x- 
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Fig. 21. — Photodissociation probability of H2 for a model («h = 100 cm -3 , \ = 1> = 1) 
computed within the FGK approximation and the exact radiative transfer in the H2 lines up 
to J = 5. 
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Fig. 22. — Photodissociation and photoionisation probability of CO and C for a model 
(rin = 100 cm -3 , x = 1, Ay = 1) within the FGK approximation and the exact radiative 
transfer in the H 2 lines up to J = 5. 
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Fig. 23. — Temperature of the gas at the edge of the PDR as a function of nu and x- Values 
next to isocontours are temperatures in Kelvin. 



-78- 




o 

-I— ' 

o 

_CD 

CD 

O 
-I—' 

o 



10 



■22 




1 = 
% = 
% = 
% = 
% = 



= 10' 



10 

10; 

= 10* 
= 10' 







10 15 
Charge of grains 



20 



25 



Fig. 25. — Photoelectric heating (erg cm -3 s -1 ) as a function of the charge of grains for 
models with n H = 1000 cm~ 3 and various values of x- Grains considered here are the ones 
with the smallest radius. The edge of the cloud is on the right of the figure where the charge 
is maximal. 
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Fig. 26. — Intensity of the 1-0 S(l) H 2 transition in 10 6 erg s 1 cm 2 sr 1 assuming a face 
on geometry. 
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Fig. 33. — Intensity of CO (2-1) in 10 6 erg cm 2 s 1 sr 1 as a function of nu and x- 
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Fig. 35. — Ratios of antenna temperatures of CO 6-5 to 2-1 transition as a function of tin 
and x- 



